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We present a new method for the computation of Lyapunov exponents utilizing repre- 
sentations of orthogonal matrices applied to decompositions of M or MM where M is the 
tangent map. This method uses a minimal set of variables, does not require renormalization 
or reorthogonalization, can be used to efficiently compute partial Lyapunov spectra, and 
does not break down when the Lyapunov spectrum is degenerate. 
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Chaotic dynamics has been investigated in a very 
large class of systems, including astrophysical, biologi- 
cal, and chemical systems, mechanical devices, models 
of the weather, lasers, plasmas, and fluids, to mention 
a few. Lyapunov exponents provide the single most im- 
portant quantitative characterization of the exponential 
divergence of initially nearby trajectories, which is the 
hallmark of chaos. Recent applications of these expo- 
nents include the connection between chaotic dynamics 
and transport theory in statistical mechanics Jl],@] and 
galactic dynamics || . 

Several methods exist for computing Lyapunov expo- 
nents However, no single method appears to 
be optimal. For example, QR and SVD (singular value 
decomposition) methods J5|,^| require frequent renormal- 
ization (to combat exponential growth of the separation 
vector between the fiducial and nearby trajectories) and 
reorthogonalization (to overcome the exponential col- 
lapse of initially orthogonal separation vectors onto the 
direction of maximal growth). The existing continuous 
versions of the QR and SVD methods also suffer from the 
additional disadvantage of being unable to compute the 
partial Lyapunov spectrum using a fewer number of equa- 
tions/operations than required for the computation of the 
full spectrum |J. Further, the continuous SVD method 
breaks down when computing degenerate Lyapunov spec- 
tra |J. The symplectic method is applicable only to 
Hamiltonian systems (and a few generalizations thereof) 
and has proven difficult to extend to systems of moderate 
size, though this is possible in principle |J. It also does 
not permit easy evaluation of partial Lyapunov spectra. 

The widespread perception that some form of explicit 
rescaling and reorthogonalization is necessary lies at the 
heart of most methods for computing Lyapunov expo- 
nents. In this Letter, we propose a general method 
which analytically obviates the need for rescaling and 
reorthogonalization. Our new method also does away 
with the other shortcomings listed above: A partial Lya- 
punov spectrum can be computed using a fewer number 
of equations as compared to the computation of the full 



spectrum, there is no difficulty in evaluating degenerate 
Lyapunov spectra, the equations are straightforward to 
generalize to higher dimensions, and the method uses the 
minimal set of dynamical variables. Since our method is 
based on exact differential equations for the Lyapunov 
exponents, global invarianccs of the Lyapunov spectrum 
can be preserved. 

The key feature of our approach is the use of explicit 
group theoretical representations of orthogonal matrices. 
This results in a set of coupled ordinary differential equa- 
tions for the Lyapunov exponents along with the various 
angles parametrising the orthogonal matrices. The sys- 
tem of differential equations is treated as an initial value 
problem and solved numerically to obtain the Lyapunov 
exponents. In the preferred variant of our method, the 
equations are only partially coupled leading to easy eval- 
uation of the incomplete Lyapunov spectrum. An in- 
teresting consequence of our methodology is the natural 
separation between "slow" (the exponents) and "fast" 
(the angles) pieces in the evolution equations. (This fact 
can be used to provide speed-up in numerical implemen- 
tations.) Since the structure of the coupled differential 
equations is of a special form, they may also turn out 
to be useful for analytic studies of evolution in tangent 
space. 

To begin, we consider an n dimensional continuous- 
time dynamical system, 



(1) 



where z = (z\, ■ • ■ , z n ) and F is a n-dimensional vector 
field. Let Z(i) = z(t) — zo(t) denote deviations from the 
fiducial trajectory Zo(t). Linearizing Eq. (Q) around this 
trajectory, we obtain 



^ = DP(zo(t),t)-Z 



(2) 



where DF denotes the n x n Jacobian matrix. 

Integrating the linearized equations along the fiducial 
trajectory yields the tangent map M (zq(<), t) which takes 
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the initial variables Z m into the time-evolved variables 
Z(t) = M(t)Z m (the dependence of M on the fiducial 
trajectory zg(f) is understood). Let A be an n x n matrix 
given by A = lim t _ >00 (MM) 1 / 2 *, where M denotes the 
matrix transpose of M. The Lyapunov exponents then 
equal the logarithm of the eigenvalues of A B . 

It is clear that M is of central importance in the evalu- 
ation of Lyapunov exponents. Its evolution equation can 
be easily derived: 

Instead of a brute force attack, our purpose is now to 
write M (or some combination thereof) in such a way 
that the resulting evolution equations are intrinsically 
well-behaved. One way to do this is to follow the ap- 
proach of Ref. and introduce the matrix A = MM. 
The evolution equation for A follows from <Q): 



dA 
dt 



DF A + A DF 



(4) 



The matrix A is symmetric and positive-definite [Q. 
Hence it can be written as an exponential of a symmetric 
matrix B fiof : A = e B . Furthermore, any symmetric ma- 
trix can be diagonalized using orthogonal matrices [fl0|| . 
Thus, A = e , where O is an n x n orthogonal 

matrix, D is an n x n diagonal matrix and O^ 1 = O. 
From standard properties of matrix exponentials, it fol- 
lows that A = Oe D 0~ 1 . There is no need for rescaling 
since the diagonal matrix D is already in the exponent 
(the diagonal elements are just the Lyapunov exponents 
multiplied by time). 

To proceed further, we use an easy to obtain explicit 
representation of the orthogonal matrix O from group 
representation theory [fi"l| . One advantage is that a min- 
imum number of variables is used to characterize the sys- 
tem: n(n — 1)/2 in O and a further n variables in D, for a 
total of n(n + l)/2. Another advantage is that numerical 
errors can never lead to loss of orthogonality. Finally, the 
dynamical equations (|i|) are solved numerically. 

We describe in detail below a variant of this idea which 
has certain further advantages. As is well-known pQj , the 
matrix M can be written as the product M = QR of an 
orthogonal n x n matrix Q and an upper-triangular n x n 
matrix R with positive diagonal entries. Substituting this 
into Eq. M), we obtain: 



QR + QR = DF QR , 



(5) 



where the overdot denotes a time derivative. Multiplying 
the above equation by Q from the left and from the 
right, we get 



QQ + RR- 1 = Q DF Q 



(6) 



Note that QQ is a skew(anti)-symmetric matrix for 
any orthogonal matrix Q and RR -1 is still an upper- 
triangular matrix. 



As before, we now employ an explicit representation 
of the orthogonal matrix Q representing it as a product 
of n(n — l)/2 orthogonal matrices, each of which corre- 
sponds to a simple rotation in the i — jth plane (i < j). 
Denoting the matrix corresponding to this rotation by 
Ow) , its matrix elements are given by: 

0%p = 1 if k = l + i,j ; 

= cos (f) if k = I = i or j ; 

= sin <j> if k = i, I = j ; 

= — sin <j) if k = j, I = i ; 

= otherwise. (7) 

Here 4> denotes an angle variable. Thus, the nxn matrix 
Q is represented by: 



Q = (9 (12) (9 (13) • ■ ■ {ln} {23) ■ ■ ■ C^™- 1 '") 



(8) 



Hence Q is parametrized by n(n — l)/2 angles which we 
denote by 9i (i = 1, • • • ,n(n — l)/2). These angles will 
be collectively denoted by 9. 

Since the upper-triangular matrix R has positive diag- 
onal entries, it can be represented as follows: 



R = 



3 Ai 



o - A 



r u ... 

e Aa r 23 



\ 



e A„ J 



(9) 



The quantities A; will be shown to be intimately related 
to the Lyapunov exponents. Our final equations will be 
in terms of the A; which already appear in the exponent, 
thus removing the need for rescaling. The quantities 
represent the supra-diagonal terms in R. 

Using the above representations of Q and R, we obtain: 



QQ = 



( -h{9) 

MO) o 

V/«-l(0) ••■ 



/n(n-l)/2(0) 



-fn-10) \ 
-/2n-3(0) 

J 



(10) 



and 



( r 'i2 r 'ln \ 



RR- 1 



A 2 r' 23 ■■■ 



' 2n 



(11) 



V A„ / 



Here, each of the n(n — 1)/2 functions fa depend (in prin- 
ciple) on the time derivatives 6\ of all the angles used to 
represent Q. In fact, they actually depend only on a 
subset of the angles. The quantities r[j are of no concern 
since they are not present in the final equations. 
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Substituting the above two expressions in Eq. (g) we 
obtain 



/ Ai ri' a r'i \ 

h{6) A 2 r 2 ' 3 • .. r 2 '„ 

\ fn-l{9) /n(n-l)/a(0) A n / 



= Q DF Q 



(12) 



Denoting the matrix Q DF Q by 5* and comparing diag- 
onal elements on both sides of ( |l2|) one gets: 



A?" — Sii 



i = l,2, 



(13) 



It can be shown || that the Lyapunov exponents are 
equal to Xi/t in the limit t — > oo. Thus, the Lyapunov 
exponents can be obtained by solving the above differ- 
ential equations for long times. However, since the right 
hand side depends on the angles 6i , we also require differ- 
ential equations governing the evolution of these angles. 

Differential equations for the angles can be obtained by 
comparing the sub-diagonal elements in Eq. ([l2]). This 
gives 

/l(0) — S21; f%{0) = S31; /n(ra-l)/2(#) = Sn,n-l- 

This set of differential equations can be transformed into 
a more convenient form Ifl2l 



gi {9), * = l,2,...,n(n-l)/2 



(14) 



where the equations for 9i are decoupled from the equa- 
tions for This avoids potential problems with degen- 
erate Lyapunov spectra. Because of these reasons, the 
second method just described is to be preferred over the 
method first discussed. Equations (|l^) and (14) form a 
system of n(n + l)/2 ordinary differential equations that 
can be solved to obtain the Lyapunov exponents. 

Our system of differential equations has another at- 
tractive feature. The equation for Ai depends only on the 
first (n — 1) 6i's (under a suitable ordering) fj^] . There- 
fore, if one is interested in only the largest Lyapunov 
exponent, one needs to solve only n equations (as op- 
posed to n(n+ l)/2 for the full spectrum). The equation 
for A2 depends only on the first 2n — 3 0i's. Therefore, to 
obtain the first 2 Lyapunov exponents, one needs to solve 
only 2n — 1 equations. In general, to solve for the first m 
Lyapunov exponents, one has to solve m(2n — m+ l)/2 
equations which is always less than n(n+ l)/2 (the total 
number of equations) for m < n. This is in contrast to 
the situation for the conventional continuous QR or SVD 
methods where it is computationally costlier to evaluate 
a partial spectrum once a threshold is crossed Q. The 
first method discussed above shares this disadvantage. 

We end the general analysis of our system of equations 
by pointing out an interesting fact. From Eq. ((L3|), 



Ai + A 2 H h A„ = Trace(S') 



(15) 



Parametrizing the Jacobian matrix DF as [DF]y = dfjj 
we can evaluate the trace of the matrix S to obtain fl2[| 

Ai + A 2 + • • • + A n = df u + df 2 2 + ■ ■ ■ + df nn . (16) 

This relation can be used to speed up numerical integra- 
tion of the differential equation for A„ . 

We now illustrate the second method for a simple two- 
dimensional example. In this case, Q is parametrized as 
follows: 



Q 



cos Oi sm tli 
— sin Q\ cos Q\ 



(17) 



and the upper-triangular matrix R may be written as, 



R = 



e Al ri2 



(18) 



e A2 

Next, we parametrize the Jacobian matrix DF as follows: 
dfu df i2 



DF 



df 21 df 2 



(19) 



Substituting the above into Eq. (|12|), we obtain the de- 
sired equations for Ai, A 2 and 6±: 

—± = df u cos 2 0i + df 22 sin 2 9 1 
dt 



-~{df 12 +df 21) sm26 x 



dX 2 



^ = d/n sin 2 0i + df 22 cos 2 0i 
dt 

+ \{df 12 +df 21) sm29 x , 
d9 1 

^ = -2( rf /n-4f22)sin20i 



(20) 



+d/i2 sin 2 0i - d/21 cos 2 0i . 

The above differential equations are numerically inte- 
grated forward in time until the desired convergence for 
the exponents, Ai/i and X 2 /t, is achieved. 

As our first example, we consider the driven van der 
Pol oscillator: 



z\ = z 2 



(21) 



z 2 = —d(l — zl)z 2 — Z\ + fccoswt 



For the parameter values d — —5, 6 = 5, and u — 2.466, 
our results are in agreement with values obtained earlier 
using the symplectic approach Q . 

To illustrate the application of the method to a system 
with more degrees of freedom, we turn to the standard 
test case of the Lorenz equations 0] : 



z\ = o(z 2 - Zl) , 
z 2 = Zx(p - z 3 ) - z 2 
h = z\z 2 - (3z 3 . 



(22) 
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For this three degree of freedom system, we need to gen- 
eralize the equations given in Eq. (|2(]). This can be 
easily done to obtain six partially coupled differential 
equations governing the evolution of the three Lyapunov 
exponents and three angles. We used parameter values 
of <t = 10, p = 28, and f3 = 8/3. An extensive com- 
parison of our method against the standard QR method 
with Gramm-Schmidt reorthogonalization (QR/GS) || 
was carried out. Both methods were applied to the same 
fiducial trajectory generated using a RK4 integrator ap- 
plied to Eqs. (23) with time step e = 0.001. Error and 
convergence analysis was carried out by applying the two 
methods to the fiducial trajectory sampled over time in- 
tervals t s > e. Both methods were implemented using 
RK4 integrators, and with t s = e = 0.001, both gen- 
erated essentially identical results. As a function of t s , 
both methods were quartically convergent as expected, 
QR/GS possessing a smaller coefficient for the positive 
Lyapunov exponent and a larger one for the negative ex- 
ponent. Even for this small system, execution times for 
both methods were similar. (We did not attempt to fully 
optimize either one of the codes.) For larger systems our 
method is expected to be more efficient. 

In the Lorenz system of equations, the sum of the three 
Lyapunov exponents must equal — (<j + j3 + 1). With 
our method, the sum of the three Lyapunov exponents 
— (cr + [3 + 1) = —13.6666 • • • was maintained to nine dec- 
imal places, independent of the sampling interval over 
the investigated range, 0.001 < t s < 0.02, a property 
not shared by QR/GS. (The sum of all Lyapunov expo- 
nents is an important quantity in stationary, thermostat- 
ted nonequilibrium systems since it is directly propor- 
tional to the transport coefficients. Recent analytic and 
numerical results are reported in Refs. Q|.) 

To summarize, we have described a technique for com- 
puting Lyapunov exponents that has several advantages 
over existing methods. The minimal number of variables 
is used, rescaling and reorthogonalization are eliminated, 
partial Lyapunov spectra can be calculated using a fewer 
number of equations, there are no difficulties with de- 
generate Lyapunov spectra, and global invariances of the 
Lyapunov spectrum can be explicitly preserved. The 
method allows a natural fast /slow split between variables 
which may be taken advantage of to improve convergence 
of the exponents. Moreover, the simple form of the fi- 
nal set of equations may prove to be useful in analytic 
considerations as well. Further details will be presented 
elsewhere |0. 
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